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1. Introduction 

Work on population ecology carried out in the 1970's certainly helped chaos to move 
center stage as a new interdisciplinary subject |1]. One- dimensional discrete-time models 
are technically among the simplest ones to consider and interpret. They provide an 
appropriate description of species with non-overlapping generations. The basic example 
for the evolution of a population density X„ G R at generation n is the linear law 
Xn+i = rXn- Here r represents the growth rate or fecundity, assumed constant. It is, 
however, too schematic allowing only extinction (r < 1), equilibrium (r = 1) or infinite 
growth (r > 1). It was soon recognized that more realistic models should be nonlinear: 



Written in this form, F is the dimensionless non-constant fitness function of the 
population or per-capita growth rate. A sage choice for it should capture the essential 
features of the system. The crucial point is that once nonlinearity is let in, a huge 
variety of new phenomena may appear as is now universally recognized. Thus, discrete 
time population ecology is in close contact with discrete dynamical system theory. For 
instance, the choice [21 [3] 



renders easy the numerical determination of the parameter values from experimental 
population data by a linear fit to logX^+i versus logX„, which constitutes certainly 
a salient advantage. In ([21) the presence of parameter K, the conventional carrying 
capacity, ensures the dimensionless character of F. A slight variant of it reads [31, 



where C is a threshold population density and the fitness parameter 6 > is 
dimensionless. 

For the purposes of the present paper we emphasize that these and similar models 
can also be used as mathematical instances of dynamical systems to illustrate different 
features when the ranges of parameters and time variable are enlarged beyond those 
realistic in population dynamics. In this spirit, we analyze here the one-dimensional 
discrete model associated with equation ^ which will be referred to as VGH (Varley- 
Gradwell-Hassell) [5]. In [6] it is briefly explained that the system is chaotic for b > 2 
and regular for b < 2, pointing out that in the order-to-chaos transition no cascade of 
period doubling emerges. In the present paper we analyze numerically and analytically 
how such a transition takes place which, to the best of our knowledge, has not been 
studied in detail in the mathematics or ecology literature. 

The system function is continuous for C = K and discontinuous otherwise. 
Piecewise continuous maps are, by far, less analyzed than continuous ones in the 
literature. Also, we address our attention to the role of the threshold parameter C, 
to which no much attention seems to have been given. 




(1) 




(2) 




(3) 
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The paper is organized as follows. In Section [2] we describe the features of the VGH 
map, present some traits of the associated dynamical system, and briefly introduce the 
pertinent equations to compute the Lyapunov exponent. In Section [3] the transition 
order-to-chaos is studied in detail. Some exact results are presented for 6 = 2 and 
an explanation for the observed bifurcation is propounded. In Section H] the effect of 
varying C and r is analyzed. Finally, Section [5] contains our conclusions. For the sake 



of completeness, two technical Appendices are included. In particular, [Appendix B 



contains a thorough analytical study of the VGH map in the critical case 6 = 2. 

2. Description of the map and its dynamics 

2.1. Alternative formulations 

We write the VGH map, equation ([3]), in the form 

Xn+i = fixn;b,c,r) (4) 

with 

f{x;b,c,r) = \ ^^[_^ ^ ^ ^ (5) 
I rx , X > c 

In going from ([1]) and ([3]) to (jl]) and ([5]) we have used the carrying capacity K as 
our natural unit to measure population densities and correspondingly introduced the 
dimensionless variable x„ = X^/K, and parameter c = C/K. 

The function / in ([5]) is defined in although, being a population variable, we 
take X G [0, oo). In the three dimensional parameter space [b, c, r) we consider only the 
region 6 > 1 to ensure decrease of / with increasing x; c,r > for compatibility with 
the range for x. 

According to this map the population density in a given generation is a linear 
function of that in the previous one as far as it does not exceed a critical value c. For 
greater values the population follows a nonlinear power-law. Figure [1] illustrates the 
different possibilities and we expect, therefore, differences in the response of the system 
according to whether the value of c is below or above unity. 

For some purposes, we have found it useful to express the VGH map in terms of 
the new variable z and the new parameter ^ 

z = 2 log(x)/ log(r), ^ = 2 log(c)/ log(r) (6) 

with X and r > 1. The numeric factor two in ([6]) has been introduced for 
convenience. The inverse transformation reads 

X = r^/\ c = r«/2. (7) 
This leads from ([5]) to 

~l" 2, ^ ^ /q\ 

{l-b)z^ + 2, Zn>i ^ ' 



z. 



n+1 
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now with phase space (—00, 00). This representation of the map has just two effective 
parameters. From a mathematical point of view this transformed version is piecewise 
hnear, whereas ([5]) defines a piecewise continuous nonhnear system. Linearization is a 
standard procedure in the study of dynamical systems. It usually follows from first order 
approximations. Here, however, the linearization is exact. In Appendix A we collect 
some results obtained from this version of the dynamical system. It is worth mentioning 
that dynamical systems of the same type can be found in applications in electronics, 
robotics and mechanical systems with impacts [71 [8] • 



2.2. VGH dynamical system 

We present here some of the most apparent features of the dynamical system originated 
by the iteration of the VGH map. 

First we can dispose of the case r < 1. a; = is then always an stable fixed point. 
When c < r^/^ there is a second fixed point at x = r^l^ which is stable if 1 < 6 < 2 and 
which attracts any initial condition in the interval (c, x*) with x* = {cjr')^!'^^^^^ . For 
6 = 2, that interval consists of period-2 points, with the exception of the fixed point 
\fr. Values x <c and a; > x* go to x = 0. Observe that here r < 1 does not necessarily 
imply extinction. This is due to the interplay between the different parameters. 

For r = 1 any initial condition in [0, c] remains fixed. If c < 1 then x = 1 is also 
fixed. It is stable for 1 < 6 < 2 attracting any initial condition in the interval (c, x*) 
with X* = c^/'-^^''^ If 6 = 2, the interior points of that interval, except x = 1, are paired 
in 2-cycles. Points to the right of x* are eventually fixed, reaching its limiting fixed 
value in two steps. 

For r > 1, X = is an unstable fixed point, independently of c. If c < r^/*, also 
X = r^/* is fixed, although it is stable only for 1 < 6 < 2. 

In the rest of the paper we consider only r > 1 which, as we shall see, generates a 
much richer behavior. 



2.3. Lyapunov exponent 

The Lyapunov exponent is a measure of the rate at which two initially close trajectories 
move away. For the one- dimensional discrete map Xn+i = f{xn) the computation of the 
Lyapunov exponent A admits an analytical formula (see for instance P [THl [HI fT2] ) 



A= hm \-J2^n\f{x,)\} (9) 
L fc=o 

which, in the case 05]), leads to 

A = lnr+ lim < - V[ln|l-6| - bln\xk\]e{xk - c) I, (10) 



fc=0 
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whereas ([H]) yields 

A = ln|l-6| lim J 1 |j ^(z, - ^ i . (11) 



k=0 

Both formulas are given in terms of Heaviside 6 function which selects the iterations 
that visit Xn > c or Zn > ^ respectively. 

The mathematical equivalence of the two expressions for A in ffTOj) and ffTTj) yields 
the following result for the statistics of points Xfc > c in the attractor 

' ^ n— 1 

-y^ln{xk)9{xk - c 



lim 

n^oo 



n 

k=0 



Inr , , 

- (12) 



which tells that the geometric mean of all the points > c in the trajectory equals the 
value r^^''. A simpler version of (ITOl) reads now 

X = \n\l-b\ Mm l-y"e(xk-c)\ . (13) 

n-*OD n ^-^ 

\ fc=0 ) 

Equation ([9]) makes use of the derivative /' which in the present case is not defined 
at x = c or 2; = ^. As it is an isolated point, it does not amount any numerical difficulty. 

In practice, one has to take care of discarding transients in order to allow the orbit 
to enter the attractor. In particular, transients are very long for chaotic trajectories 
when 6^2. Besides, it is convenient to average A over a number of initial seeds. 

The interpretation of the Lyapunov exponent in this system is particularly 
straightforward: A is ruled by the proportion of exterior points in the trajectories, 
namely those with x„ > c or > ^. 

Since the summation in (ITOl) and (ITTl) is over a large enough number of points 
on the trajectory, numerical accuracy is a relevant issue. To this end, and following 
the analysis in [13] , we have used the multi-precision Fortran package MPFUN p3l [T5] 
which allows one to efficiently compute with very high number of digits. We concluded 
that double precision gives here enough computational guaranties as regards (fTOl) and 
(ITTl) . By contrast, it is not the case for the computation of trajectories of (IHl) with some 
integer values 6 > 2, an issue that we deal with in Section 14.51 



3. Dependence on 6: order-to-chaos transition 

As already stated in |6j, the chaotic regime corresponds to 6 > 2, independently of the 
values of c and r, without stable regularity windows embedded. This is clear from the 
positive character of A in f|T3l) when b > 2. 

In the limit 6-^2+ the Lyapunov exponent becomes very small due to the 
term ln|l — b\, since the proportion of exterior points is finite. Following results on 
nonextensive thermodynamics it may be pertinent to wonder whether a power law n'^ 
is a more appropriate ansatz than the exponential exp(A?T,) to measure the divergence 
of initially close trajectories fiQ[ ITTj . We have carried out a statistics of the largest 
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separations, in this limit, as a function of n and the results point out towards a clear 
exponential divergence. 

The behaviour of the system with respect to b is illustrated by the bifurcation 
diagrams in Figure [2j The three panels correspond to c = 1.2, 1, 0.8, from top to 
bottom, respectively, with the fixed value r = 2. In Figure [3] we plot the bifurcation 
diagram for c = 2.5 and r = 2, as well as the Lyapunov exponent A (upper panel). 
Notice the irregularities in the curve around the value 6 = 2.2 where the bifurcation 
diagram exhibits a number of band merging crises (121 [18] . 

The aforementioned sudden transition from a regular to a chaotic system, with no 
period doubling, is apparent. As the negative sign of the computed Lyapunov exponent 
witnesses, b < 2 means regular behavior: every initial seed ends up in a periodic orbit. 

For fixed b < 2, the attractor consists of coexisting limit cycles. Its cardinal is a 
piecewise constant function of c. So, there exist some values c = Cn {n > 1 integer) 
that punctuate the sequence of the number of points that make up the attractor. The 
cardinal is given by the heuristic formula 2{n + l)/[3 — (—1)"], (n > 1), whose first few 
terms read 

1, 3, 2, 5, 3, 7, 4, 9, 5, 11, 6, 13, 7, 15, 8, 17, 9, . . . (14) 

They can be interpreted as forming two intermingled sequences whose terms increase 
by one and two units respectively. 

The diagram in Figure H] gives a perspective of this situation for b = 1.9 and r = 2, 
as a function of c. 

In the rest of this Section we focus our attention on the case 6 = 2 which separates 
the regular and chaotic regimes. We present both numerical and analytical results. The 



reader is referred to Appendix B for details. 



3.1. Case 6 = 2. 

For 6 = 2 and fixed r > 1 the orbit of an arbitrary initial condition xq G [0, oo), is simple 
enough to be rigorously described. The behavior is varied and shows an interesting 
dependence on the parameter c. In particular this permits to appraise the effect of the 
discontinuity. As detailed below, the general characteristic in this case is that all the 
initial conditions are fixed, periodic or eventually periodic. There exist an infinity of 
coexisting limit cycles. Their periods are allowed to take only a restricted set of values, 
which depend on the initial condition and parameters value through very strict laws. 



We defer their detailed description to Appendix B The reader can find there precise 



information about the transients length and exact account of the cycles. Here we report 
only the most salient features in terms of the variable x and threshold parameter c. 

If c < 1, then the closed interval [c, r/c] is invariant under the action of the map. 
Any interior point xq ^ ^/r belongs to the 2— cycle {xq, t/xq}. Xq = -y/r is a fixed point. 
The extremes Xq = c and Xq = r/c as well as all the exterior points are eventually 
periodic, entering the invariant interval after a transient. 



Dynamics of a map with power-law tail 



7 



If c > 1, then the closed interval [1/c, rc] is invariant under the action of the map. 
Every exterior point is eventually periodic. Interior points are periodic. For fixed c, the 
period of the cycles depends on xq but can take on only a restricted set of values. For 
instance, with c = 2 one finds a 3— cycle: {1/y/r, ^/r,r^/r}■, a 4— cycle: {1/r, 1, r, r^}; 
and the rest of the points in the interval accommodates in an infinity of 6— cycles. All 
these cycles are the limit cycles for the eventually periodic points. Notice that in this 
system the existence of period three does not imply chaos. This is not in conflict with 
the celebrated Li and Yorke theorem because the VGH map is defined by a discontinuous 
function. 

The basins of attraction of such limit cycles are infinite intermingled sets of zero 
measure made of equidistant points in the z variable. The vertical segments located 
at 6 = 2 in Figures [2] and [3] comprise all the coexisting limit cycles. In numerical 
simulations such a vertical segment gets filled only if enough initial conditions are used 
in the construction of the bifurcation diagram. 

3.2. An explanation of the observed bifurcation 

Next we develop an heuristic explanation for the origin of the order-to-chaos bifurcation 
observed at 6 = 2. As is well known, period-2 points are fixed points of the second 
iterate /^^l of the map. In the cobweb plot of /I^' shown in Figure [5] the relevant element 
is the piece of /'^l which is co-linear with the bisectrix, because it conveys the existence 
of a continuous interval of fixed points. Leaving apart the fixed point of /, these points 
pair themselves to form the alluded period-2 orbits. By the cobweb construction the 
rest of the points eventually land in the piece on the bisectrix. 

Whenever an infinite set of trajectories of period n > 2 is observed, with 6 = 2, it 
is because this value of the parameter forces /I""' to have one or more segments co-linear 
with the bisectrix too. The question is then: Which are the conditions for the map / 
to allow iterates /'"^ to have pieces co-linear with the bisectrix? 

Our answer is that the phenomenon may take place whenever the piecewise 
continuous map / has a piece defined by a power-law function. To buttress this answer 
let us write the general expression of a term of /f"' as 

rPib)^ii~br (15) 

where p{b) stands for a polynomial of b and m < n E N. Since 6 > 1, the conditions for 
( |T5l) to be co-linear with the bisectrix are: even m, b = 2 and p{2) = 0. 

To further illustrate this, let us focus on one part of the second iterate of ([5]), namely 



where the dots stand for the three remaining terms and intervals in the definition of 
the second iterate. The explicit term of /[^^ written down in f|T6|) is co-linear with the 
bisectrix in the cobweb plot if and only if 6 = 2, as expected, provided c < ^/r. If 




(16) 
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c > -y/r the extremes of the interval collapse to a single point and that particular piece 
disappears in the second iterate. 

In principle, several iterates can have co-linear terms simultaneously for the same 
value of c. This gives rise to coexisting cycles of different periods, a fact that is illustrated 
in the bifurcation diagram of Figure El There we have coded in color (gray tone) 
the periods of the trajectories, with fixed values b = r = 2. The representation is 
given in terms of {z,C,}, which endows the plot with more symmetric structure. The 
intermingling of the colored rhombuses accounts for the coexistence of limit cycles with 
various periods. In addition, horizontal dashed (color) segments stand for cycles whose 
periods are indicated by circled numbers and have been taken from Table IBli Framed 
numbers indicate the corresponding period of orbits in rhombuses. As Figure [6] shows, 
for the parameter values considered, only two different iterates can have co-linear terms 
for the same value of ^. Moreover, as ^ increases, more co-linear terms of successively 
higher iterates emerge whereas the co-linear terms of the previous iterates collapse and 



disappear (see Appendix B) 



In other words, the presence of the rx^~^ power-law piece in the definition of the 
VGH model is at the origin of the bifurcation at 6 = 2. More generally, any one- 
dimensional map defined as a piecewise function, with a power-law piece in its definition, 
is a candidate to exhibit a bifurcation of the present type. 

3.3. Case 6 = 2. When xq = c: Harter's boundaries 

As is well known, for a continuous unimodal uniparametric map with its maximum 
at X = X*, the plot of the successive iterates of the initial seed xq = x* versus the 
map parameter generates the so-called Harter boundaries [HI [20] . In the logistic map 
[H [m [211 [22] , instance, they correspond to the sharp cusps observed in the invariant 
density. These boundaries are, to some extent, the skeletal system of the bifurcation 
diagram where they are readily observed provided it is built up using color or grey tones. 
Harter lines cross at unstable equilibrium points and, in the case of the logistic map, 
the set of crossing points follows itself a reversed bifurcation cascade. 

The VGH map is not differentiable at x = c, although it is still unimodal in the 
sense of being monotonically increasing for x < c and decreasing for x > c. A similar 
study to the one described for the logistic map leads to interesting results. To be precise, 
Harter lines correspond to the color boundaries visible in Figures [2 and [3l They can 
be directly established from the form of the map. It is nevertheless simpler to use 



instead ([8]) and the results gathered in Appendix A For b = 2, the seed xq = c ends up 



in a cycle whose period P(r, c) depends on r and c. The expression for P and the cycle 
elements may be explicitly written down. 

For c < 1 the 2-cycle {rc, 1/c} is reached just after one iteration. For c > 1, x = c 
is always a periodic point. Obviously c = 1 belongs to the 2-cycle {1,t}. When c > 1 
there is a unique integer M G N such that 
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with 



M 



jlogc 
logr 



in terms of the integer part or floor function. Then, for c > 1 we get two cases 



P(r, c) = M + 2, if c = r*^/2 



P{r, c) = 2(M + 2), if c G {r^'^^ r'^^^+m^^ 
The 2(M + 2)-points of the cycle for the case (l20l) read 

C C C J. J.2 ^M+l 



which in turn, for c 



r 



Af/2 



[i.e., case (1191) ) contract to the (M + 2)-points cycle 



{ 



-M/2 



,r 



M/2-1 ^Af/2+1 



r 



} 



(18) 

(19) 
(20) 

(21) 

(22) 



It is interesting to observe that fl2Tl) can be obtained recurrently from {rc, 1/c} by 
a simple procedure. In going from M to M + 1 two new elements are added to the 
cycle: the one in the leftmost position is obtained by dividing by r the first element 
in the previous cycle. The other, which goes at the rightmost position, is obtained by 
multiplying by r the last element of the previous cycle. 



4. Dependence on c: Discontinuity location 

The dependence of the system with respect to the parameter c with r fixed, or 
equivalently ^, exhibits a large variety of features. We describe the system in the chaotic 
regime. First near the bifurcation point 6=2, and next for higher h. 

4.1. Nearb = 2 

In Figure [7] we plot the bifurcation diagram as a function of c, computed with b = 2.01 
and r = 4. For c < 1 the trajectories Xn wander in two relatively large chaotic bands in 
contrast with the extremely narrow ones for c > 1. 

The upper panel in the figure gives the Lyapunov exponent. It is everywhere 
positive, as it must be for a chaotic regime. As a function of c, its value changes 
suddenly at every band merging crisis observed in the bifurcation diagram. 

To study the diagram with further detail we present in Figures IHl-fTOlmagnifications 
of small areas in Figure [3 

The region with c < 1, near 6 = 2, has further interesting features. Thus, what 
in Figure [7] appears as two dense chaotic bands does have structure when minutely 
examined. Figure M shows this feature for the upper band. The grid structure in the 
bifurcation diagrams as a function of c (or C,) fades as c decreases. The Lyapunov 
exponent in the upper panels exhibits jumps associated with chaotic bands merging 
crises. 

For c > 1, the bifurcation diagram in Figure [8] exhibits further details after 
magnification. A zoom of the area indicated by the arrow is given in Figure [9] where a 
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sudden variation of the size of the attr actor is apparent. This thin window is composed 
by twenty two chaotic narrow bands (only five in the zoom). Further, the crossing of 
Harter hues in Figure [TT] punctuate the start and end of the shrunk chaotic bands. A 
unstable orbit exists inside the window, which is represented by the dashed lines in the 
Figure [Til Consequently, what we observe in Figure [9] as a shrinking and widening of 
the attractor corresponds, actually, to a pair of interior crisis [121 UHl 123]. Moreover, 
further interior crises of various smaller sizes appear inside the window itself in a, 
possibly, self-similar way. For instance, in the leftmost part of this bifurcation diagram 
(1.14 < c < 1.1405) one can hint a crisis of the very same kind. 

Figure [THl corresponds to a zoom of the left and uppermost part of the bifurcation 
diagram in Figure [81 The Lyapunov exponent in the upper panel shows oscillations at 
band merging crises and, once again, a plateau (only visible in the inset) at the crisis 
located at 0.963 < c < 0.964. An explanation for the occurrence of A plateaus is deferred 
to Section [T3l 

For c > 1 and close to the critical point 6 = 2, the crises take place at integer values 
of ^. This feature is illustrated in Figure [T2l The left panel shows A as a function of c for 
b = 2.01 and three different values of r. The right panel corresponds to the same data 
expressed in terms of ^: All three curves collapse onto a master curve. Furthermore, 
the value of the Lyapunov exponent is invariant also in magnitude. This property fades 
as far as we move toward higher values of b, namely far from the critical point. All this 
buttresses the existence of universality in the system near the transition order-to-chaos. 

Bifurcation diagrams of the VGH system always collapse when expressed in terms 
of z versus ^. This is true even far from 6 = 2, in contrast with the Lyapunov exponent 
diagrams. 

4 ■ 2. Far from 6 = 2 

The lower panel of Figure [131 gives a view of the bifurcation diagram near c = 1, with 
r = 4 and b = 2.1, i.e. leaving the critical point. The upper panel shows an enlargement 
of the region located by the arrow below. This is again an instance of interior crisis. 
Concomitantly, the uppermost panel shows the variation of the Lyapunov exponent. 
At variance with Figures [SHTOl here A exhibits neither jumps nor plateaus across the 
window where the attractor shrinks. The reason is explained in the next subsection. 

In the bifurcation diagram of Figure with b = 2.2 and r = 2, the band merging 
phenomenon is more involved than in Figure [3 The corresponding Lyapunov exponent, 
in the upper panel, exhibits large variations at crises too. However, these do not take 
place any longer at integer values of ^. 

Eventually, in Figure [151 where 6 = 3, we observe that the Lyapunov exponent 
presents a large plateau whereas in the corresponding bifurcation diagram no interior 
crisis occurs. We have not found a justification for this case. 
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4-3. Interior crises and X plateaus 

An plausible explanation for the origin of A plateaus occurring between interior crises 
pairs reads as follows. According to our numerical simulations, the narrow windows 
between interior crises of the VGH map resemble very much the so-called cycles of 
chaotic intervals, where the trajectories jump in a cyclic way from one chaotic band to 
another. Such a phenomenon occurs, for instance, in the logistic map [21] and is also 
termed in the literature as cyclic chaotic attractor or cyclic chaotic bands |25j . 

Trajectories look quite regular and, if the line x„ = c in the bifurcation diagram 
does not cross any of the thin bands inside the window, then the proportion of exterior 
points {xn > c) remains almost constant. As a consequence of it and taking into account 
the interpretation of A for this map at the end of Section 12.31 the Lyapunov exponent 
gets the plateau shape. Otherwise, A varies across the window, which is the case in 
Figure [H 

4.4- Crisis-induced intermittency 

Next we gather some results, obtained from numerical simulations, concerning the 
behavior of trajectories around interior crises. 

We commence by pointing out that interior crises in the VGH map appear in pairs. 
As the parameter c increases through a star value cl the chaotic attractor suddenly 
shrinks into thin chaotic bands. This is the first interior crisis. Then, it exists a second 
star value > cl where the set of thin bands suddenly widen and the attractor at c > C2 
recovers its old size. This is the second interior crisis of the pair. By contrast, in the 
logistic map, a tangent bifurcation precedes always an interior crises [12]. 

For values of c slightly different than a star value, in the region where the attractor 
widens (c < cl,c > C2), the orbits spend long stretches in the region where the 
attractor is confined between the two crises. Occasionally, the trajectories burst and 
visit the whole attractor. This behavior, termed crisis-induced intermittency [121 126] , 
is illustrated in Figure [161 There, two orbits are plotted for slightly different values 
Cl = 0.8840 and C2 = 0.8845, located at both sides of the star value c* = 0.88411175 . . .. 
The empty square symbols stand for a trajectory in a regime where the attractor is 
still large. Intermittency is clearly observed. The solid dots represent a trajectory in a 
regime with shrunk attractor. It seems, at first glance, a period-8 orbit but the inset, 
which is a zoom of just one single narrow band, allows us to illustrate its non-periodic 
character. 

For fixed c, the statistics of the length of stretches where the orbit stays confined 
in the region of the shrunk attractor is well described by an exponential distribution, 
provided |c — c*| is small, which yields a characteristic length t(c). It has been shown 
[26] that for a large class of dynamical systems which exhibit crises, the dependence 
reads 



(23) 
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Indeed, this is the case for the VGH map. We have built up the statistics of r as a 
function of |c — c*| for the case in Figure [T6l The results are given in Figure [T71 The 
linear fit in log-log scales yields a determination of the critical exponent: 7 = 0.89±0.02. 

4-5. A numerical flaw 

Version (IHl) of the map presents a worth mentioning numerical nuisance. Namely, for 
some integer values of 6 > 2, computer generated trajectories collapse to an unstable 
periodic orbit after some iterations. A similar nature numerical artifact has been studied 
in [271 EH]. 

The dynamics for real z may be viewed as follows. For an initial point Zq < ^, every 
iteration conveys a shift to the right by two units until the condition Zn > ^ is fulfilled. 
The second line in (|8]) may be read as: {(1 — 6) [z„J + 2} + (1 — mod 1). Thus, for 
integer b > 2 the quantity in brackets is an integer and hence, its successive iteration 
conveys just a shift, as above. The key point is that, for finite machine precision, the last 
term dramatically looses precision in each iteration for some particular integer values of 
b>2. 

This effect may be understood on the basis of a generalization of the Bernoulli or 
binary shift map: Zn+i = 2zn (mod 1). To this end, let us use the binary representation 
of the number z^mod 1. Multiplication in base 2 is very simple in some cases. For 
instance, when \ l — b\ = 2^ we get the binary representation of the product 2^(2;^ mod 1) 
simply by shifting k — 1 places to the left every bit and (due to the finite precision of 
the computer) adding simultaneously k — 1 zeros to the right of the number. Henceforth 
this number gets shifted by two units in every iteration till it reaches Zn > C, when 
the mod operator acts again. This way, the piece in bracket is preserved as an integer 
under the action of the floor function, whereas the term 2'^(z„ mod 1) looses significant 
digits before going through the loop again. Eventually it stops when, after a number 
of iterations, 2''{zn mod 1) = 0. At this point the iteration reduces to a periodic orbit 
with elements located at integer numbers. The full output is then a set of limit cycles 
whose elements are always integers. This misleading result is just consequence of the 
computer finite precision. 

5. Discussion and conclusions 

We have revisited a one-dimensional population model, proposed as an instance of 
density dependent dynamics. The fitness parameter b controls the onset of chaos. Values 
b > 2 convey chaos. Otherwise the system is regular. On the critical point 6 = 2 the 
attractor is made of an infinity of limit cycles that share a finite number of different 
periods. Hence, the transition order-to-chaos takes place through three steps: i) Finite 
number of limit cycles at 6 < 2; ii) Infinite number of limit cycles with finite number 
of different periods at 6 = 2; and Hi) Chaos at 6 > 2. No period doubling occurs in 
this route to chaos. At 6 = 2, we have given an exact description of the system. We 
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have checked that such regular attractors have been sometimes plotted in bifurcation 
diagrams in the literature. However we have not found the concurrent descriptions of 
the attractor. Moreover, bifurcation diagrams may be found in the literature where this 
attractor lacks, most likely because too few initial seeds were used in the computation. 
An explanation for the phenomenon to occur has been provided on the basis of the 
power-law tail of the map. 

The behavior of continuous piecewise linear maps has been studied in terms of the 
so-called border collision bifurcations. Thus, a thorough classification of the different 
kinds of bifurcations, in the parameter space, is given in [29]. The study explicitly 
excludes the borderline systems where a slope equals ±1. Interestingly, the VGH map 
in its version ([8]) with ^ = (c = 1), falls just into this category. A study of discontinuous 
maps [30], namely with ^ 7^ (or c 7^ 1), is flawed by the same restriction. Hence, the 
VGH map corresponds to a family of functions that are, by definition, out of those 
analyses. Moreover, in view of the results given in Section 13.21 this type of bifurcation 
may endow that classification with further structure. 

An interesting phenomenon occurs close to the critical value, i.e. for 6 ~ 2. 
There, the Lyapunov exponent of the VGH model exhibits renormalization features 
when expressed as a function of the variable ^. In turn, the bifurcation diagrams as a 
function of c collapse for all r and b when expressed as {z^,^}. 

The abrupt variation of the Lyapunov exponent near band merging crises has 
already been reported [31]. Moreover, we have observed plateaus in the Lyapunov 
exponent, when represented as a function of c. It is worth mentioning that interior 
crises in the VGH map appear always in pairs. 

We think that the the presence of a power-law piece in maps like the VGH model 
studied here originates worth knowing properties. We do hope that the variety of results 
gathered in this work may stimulate further work. 
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Appendix A. Exact piecewise linearization 

Although by no means necessary, version ([8]) of the VGH dynamical system is 
extremely useful for analytic purposes and has a simple mathematical interpretation: it 
corresponds to two coupled afiine djTiamical systems. Precisely the new parameter ^ 
introduced acts as coupling constant. 

For the general affine discrete dynamical system 



Xn+l = AXn + B 



(A.l) 



the general solution with initial condition xq is 




(A.2) 



Dynamics of a map with power-law tail 



14 



By itself system (lA.ip has a rather dull dynamics. The evolution tends to the fixed 
point B/{1 — A) if 1^41 < 1, whereas escapes to infinity if |y4| > 1 or A = 1. For A = —1 
any initial condition xq except the fixed point B/2 belongs to the 2-cycle {xq, —xq + B}. 

For the system ([8]) , the coupling parameter ^ divides the phase space for our variable 
z in two regions: z < ^ (Region I) and ^ > ^ (Region II). In I (lA.ll) holds with A = 1, 
whereas ^4 = 1 — 6<0in region II. All over the phase space B = 2. The previous 
equation flA.2l) is instrumental in discussing the transitions to and fro between both 
regions. These transitions explain the more complicated dynamics of the coupled system 



N 



(A.3) 



Transitions I ^ II from a point Zo < C will necessarily take place after A^ iterations 

with 

2 

independently of the value of parameter b. Furthermore the orbit of zo < will land in 
II always with & + 2]. 

In contrast, if > ^ transitions II I are not always possible and, when they 
actually occur, follow a more complicated pattern depending on b, ^ and zo- Concerning 
parameter b it is convenient to distinguish two possibilities: 

B : b>2 . ^ ^ 

As ^ and zq are concerned we set apart three cases: 
a: ^ < -20 < 

b: ^ < < ;Zo (A.5) 
c: < ^ < Zq . 

We have the following variants for case A in (IA.4p : 



Aa: transitions II ^ I are forbidden. 

Ab: transitions allowed only if 2:0 > (^ — 2) /(I — 6) and then A^ = 1 

Ac: transitions allowed for any 2:0 if > 2/6 or if ,^ < 2/6 < and zq > (^ — 2)/(l — 6). 
In any case A^ = 1. 

To discuss case B in flA.4p we introduce 
ln|(6e-2)/(6;.o-2)| 

P = Hb^) • ^ ^ ^ 

We have the following variants: 

Ba: transitions are always allowed and A^ = 2[p/2], in terms of the ceiling function. 
Bb: If Zq < 2/6 then A^ = 2[p/2]. For zq > 2/6, if zq < 4/6 - ^ then A^ = [pj + 1 or 

A^ = [pj + 2 according to [pj being even or odd. If zq> 2/b and 2;o > 4/6 — ^ then 

A^ = 1. 

Be: the trajectory always crosses to region I with varying N according to the following 
scheme: 
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(a) For ^ < ;zo < 2/6, = 2\p/2] 

(b) For ^ < (2/6) < Zq < (4/6 -0, then = [pj + 1 if [pj is even, or = [pj +2 
if LpJ is odd. 

(c) For { < 2/6 < 4/6 - ^ < 2o, = 1 

(d) For 2/6 < ^, A^ = 1, independently of 2:0 > ^ 



Appendix B. Detailed study of the system 6 = 2 

Here we obtain the basins of attraction, periods, cycles and transients for the VGH map 
with 6 = 2 and arbitrary ^ and xq- 

We analyze separately the cases ^ < (c < 1), ^ = (c = 1) and ^ > (c > 1) 
with fixed 6 = 2. To facilitate the notation, in this Appendix we will use the 2;— version 
of the map, namely 

Appendix B.l. ^ < (or c < 1) 

The first general feature of the system is that the interval 2 — ^] is invariant under 
the action of the map. This is readily appreciated on the bifurcation diagrams in figures 
E] and [3l Every initial zq in the invariant interval belongs to a cycle of period two, with 
the exception Zq = 1 which is a fixed point. The point a;o = (or equivalently Zq = —00) 
is fixed Vc (respectively V.^). The remaining points zq are eventually periodic and, after 
a transient, they enter the interval 2 — The phase space get partitioned as 

(-00, 00) = (-00, |J(^' 2 - |J[2 - 00), (B.2) 
which corresponds to the original x— phase space partition 

(0,oo) = (0,c]|J(c,^)|J[^,oo). (B.3) 
The dynamics of the map encompasses three cases, defined by the value of ^: 

(i) Zq < C, (or xo < c). Define N = [{C^ — zo)/2\. Transient: A + 1. Limit cycle: 
{2(A+l) + zo,-2A-2o}. 

(ii) ^ < 2:0 < 2 — ,^ (or c < a;o < r/c). No transient. Cycle: {zo,2 — zq}. zq = 1 fixed 
point. 

(iii) 2 — ^ > (or ^/c > Xq). Define A = [(^ + -2o)2j. Transient: A + 1. Limit cycle: 
{2{N + l)-zo,-2N + zo}. 

Appendix B.2. = (or c = 1) 

The invariant interval under the action of the map is [0,2], and their points belong to 
cycles of period two, with the exception Zq = 1 which is a fixed point. Points zq ^ [0, 2] 
are eventually periodic. The partition of the 2— phase space is 

(-00, 00) = (-00, 0) |J[0, 2] |J(2, 00), (B.4) 
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which corresponds to the original x— phase space partition 

(0,oo) = (0,l)|J[l,r]U(^'°°)- (B-^) 
We distinguish four cases: 

(i) zq < (or xo < c = 1). Define N = —[2:0/2]. Transient: N. Cycle: 

{2N + zo,-2{N+l)-zo}. 

(ii) Zq = and Zq = 2 constitute the cycle: {0, 2}. 
(in) < Zq < 2 with Zq 7^ 1. Cycle: {zq, 2 — Zq}. 

(iv) Zq > 2 (or xq > r). Define N = [2:0/2]. Transient: + 1. Limit cycle: 
{2(iV+l)-^o,^o-2Ar}. 

Appendix B.3. ^ > Q (or c> I) 

This is the most involved situation. Here the invariant interval under the action of 
the map is [—^,2 + The remaining points are eventually periodic. The convenient 
partition of the phase space is 

(-00, 00) = (-00, -0 UK' 2 + i] y (2 + 00). (B.6) 

In the original x— phase space we have the partition 

(0,oo) = (0,i)U[i rc]U(rc,oo). (B.7) 
Next we describe the features of the various intervals: 

(i) Zq < (or Xq < 1/c). Transient: [— (^ + Zq)/2\ + 1. 
(n) zo>2 + ^ (or Xq > rc). Transient: [{zq - ^)/2j + 1. 

(in) —{< ^0 < 2 + { (or 1/c < Xo < rc). This case embraces five possibilities according 
to the value of ^: 

(a) < ^ < 1. The convenient partition of the 2;— phase space is 

K,2 + e]= (B.8) 

K,o)U[o,e]U(^,2-0 

U[2-e,2)U[2,2 + ^] 

The subinterval I = (^,2 — ^) is, in turn, invariant under the action of the 
map. 2:0 = 1 is a fixed point. Else, orbits are the period-2 cycles: {2:0, 2 — 2:0}. 
The other four subintervals together form also an invariant subset. Their 
points hop in period-4 trajectories following the symbolic cyclic sequence: 

(b) ^ = 1. The interval under scrutiny is now (—1,3). There exist the cycles: 
{—1, 1, 3} and {0, 2}. Points Zq ^ —1, 0, 1, 2, 3 are in cycles of period four. 
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(B.9) 



U(0,2-OU[2-e,e]U(^,2) 
|J{2}U(2,4-OU[4-^'2 + e] 



There is no invariant subinterval under the map action. Four different types 
of periodic orbits exist. One period-2 trajectory which is the cycle {0, 2}. 
One period-3 trajectory which is the cycle { — 1,1,3}. Period-4 orbits are, 
symbolically: B ^ S C ^ J-', cyclic. Similarly, period-6 orbits follow the 
pattern: 

A^T>^Q^A^T>^Q, cychc. Notice that every one of the three 
subintervals is visited twice in one cycle. 

(d) ^ = 2. One period-3 cycle: {-1,1,3}. One period-4 cycle: {-2,0,2,4}. The 
remaining orbits are period-6. 

(e) ^ > 2. This is by far the most involved situation. We define the following 
quantities: = [^/2j and a = ^ — 2N. The convenient partition now depends 
on and a: 



The length of the subintervals J's and if s is a and 2 — a respectively. The 
super- index tells whether the interval is located either to the left (— ) or to the 
right of 2; = 0. The sub-index codes the borders location of the interval 
according to 



[-e, 2 + e] = 



(B.IO) 




{2k, 2k + a), 
{2k + a,2{k+l)), 
{-2k - a, -2k), 
{-2{k + l),-2k-a). 



Besides 



C 



A 



B 



{2k + a ; k = O...N + 1}, 
{-2k -a; k = O...N}, 
{2k; k = -N ...N + 1}, 
{2k + l; k^-N...N}. 
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Table Bl. Period of the trajectories according to the initial point zq G [~£,,2 + 
with 6 = 2 and ^ > 2, in equation ([8]). 



Initial point 


T] 


Period 


Condition 


Zq t Jj^ and (fi u 


— AK ± Zq 


4A^ + 6 
AN + A 


r] + a> 2 
1] + a <2 


zo e and ^ P 


2{k + l)TZo 


AN + 2 
AN + A 


T] > a 
1] < a 


zoeA[jB 




4A^ + 6 
AN + A 


a> 1 
a <1 


zoeV 




2N + 3 
2N+1 


a> 1 
a <1 


zoeC 




2N + 2 


Va 



The period of any trajectory starting in [— ^, 2 + ^] is given in Table IBll where 
we have defined an auxiliary quantity rj when appropriate. 



May R M 1976 Simple mathematical-models with very complicated dynamics Nature 261 459-67 
May R M, Conway G R, Hassell M P and Southwood T R E 1974 Time delays, density-dependence 

and single-species oscillations J. Anim. Ecol. 43 747-70 
Hassell M P 1975 Density-dependence in single-species populations J. Anim. Ecol. 44 283-95 
Varley G C, Gradwell G R and Hassell M P 1973 Insect population ecology: an analytical approach 

(Oxford: Blackwell) 

May R M and Oster G F 1976 Bifurcations and dynamic complexity in simple ecological models 
Am. Nat. 110 573-99 

May R M 1975 Biological populations obeying difference equations: Stable points, stable cycles, 

and chaos J. theor. Biol. 51 511-24 
Banerjee S and Verfghese G C 2001 Nonlinear phenomena in power electronics: Attractors, 

bifurcations, chaos and nonlinear control (New York: IEEE Press) 
Hogan S J, Higham L and Griffin T C L 2007 Dynamics of a piecewise linear map with a gap Proc. 

R. Soc. A 463 49-65 
Hilborn R C 2000 Chaos and nonlinear dynamics (Oxford University Press) 
Strogatz S H 1994 Nonlinear Dynamics and Chaos (Cambridge MA: Perseus Books) 
Collet P and Eckmann J P 1980 Iterated Maps on the Interval as Dynamical Systems (Birkhauser) 
Ott E 2002 Chaos in dynamical systems (Cambridge University Press) 

Oteo J A and Ros J 2007 Double precision errors in the logistic map: Statistical study and 

dynamical interpretation Phys. Rev. E 76 036214-22 
Bailey D H 1993 Algorithm 719: Multiprecision translation and execution of Fortran programs 

ACM Trans. Math. Software 19 288-319 
Bailey D H 2005 High-precision floating-point arithmetic in scientific computation Comput. Sci. 

Eng. 7 54-61 

Costa U M S, Lyra M L, Plastino A R and Tsallis C 1997 Power-law sensitivity to initial conditions 
within a logisticlike family of maps: Fractality and nonextensivity Phys. Rev. E 56 245-50 

Tsallis C, Plastino A R and Zheng W M Power-law sensitivity to initial conditiions: New entropic 
representation 1997 Chaos, Solitons & Fractals 8 885-91 



Dynamics of a map with power-law tail 



19 



Grebogi C, Ott E and Yorkc J A 1983 Crises, sudden changes in chaotic attractors, and transient 

chaos Physca D 7 181-200 
Jensen R V and Myers C R 1985 Images of the critical points of nonhnear maps Phys. Rev. A 32 

1222-4 

Eidson J. Flynn S, Holm C, Weeks D and Fox R F 1986 Elementary explanation of boundary 
shading in chaotic-attractor plots for the Feigenbaum map and the circle map Phys. Rev. A 33 
2809-12 

Feigenbaum M J 1978 Quantitative universality for a class of non- linear transformations J. Stat. 

Phys. 19 25-52 

Feigenbaum M J 1979 Universal metric properties of non-linear transformations J. Stat. Phys. 21 
669-706 

Lai Y C, Grebogi C and Yorkc J A 1992 Sudden change in the size of chaotic attractors: How does 

it occur? Applied Chaos (chapter 19) Ed. by Kim H and Stringer J (John Wiley & Sons, Inc.) 
Beck C and Schlogl F 1993 Thermodinamics of chaotic systems (Cambridge University Press) 
Maistrcnko Yu L, Maistrcnko V L and Chua L O 1993 Cycles in chaotic intervals in a time-delayed 

Chua's circuit Int. J. Bifurcation and Chaos 3 1557 72 
Grebogi C, Ott E, Romeiras F and Yorke J A 1987 Critical exponents for crisis-induced 

intermittency Phys Rev E 36 5365-79 
Diamond P, Klocdcn P, Pokrovskii A and Vladimirov A 1995 Collapsing effects in numerical 

simulations of a class of chaotic dynamical systems and random mappings with a single attracting 

center Physica D 86 559-71 
Yuan G and Yorke J A 2000 Collapsing of chaos in one dimensional maps Physica D 136 18-30 
Banerjee S, Karthik M S, Yuan G and Yorkc J A 2000 Bifurcations in one-dimensional piecewise 

smoth maps: Theory and applications in switching circuits IEEE Trans. Circuits Syst.-I: Fund. 

Theory Appl. 47 389-94 

Avrutin V and Schanz M 2006 On multi-parametric bifurcations in a scalar piecewise linear map 
Nonlinearity 19 531-52 

Mehra V and Ramaswamy R 1996 Maximal Lyapunov exponent at crises Phys. Rev. E 53 3420-4 



Dynamics of a map with power-law tail 



20 




I I I I L 

1 2 

X 



Figure 1. Shape of the map (O for 6 = 2.5, r — 2 and three different values 
c — 0.7, 1, 1.3. The curves have been vertically shifted for the sake of clarity. 




Figure 2. (Color online) Bifurcation diagrams as a function of b, for r — 2 and three 
different values of c. The color scale is logarithmic and it stands for the frequency the 
point is visited with. This holds for the rest of color figures, unless otherwise stated. 
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Figure 4. (Color online) Regular orbits. Bifurcation diagram as a function of c, with 
b — 1.9 and r — 2. The period of cycles is coded by the color scale (gray tone). 



Dynamics of a map with power-law tail 



22 




Xn 



Figure 5. (Color online) Cobweb plot for /'^l (thick red line). Parameter values: 
b = r = 2, c ~ 0.8. The three seeds are Xq = 0.11,0.7,2.6, and are eventually stable 
fixed points of f^"^^ . 
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-10 1 2 3 

Figure 6. (Color online) Bifurcation diagram for r = b = 2. The color (gray tone) 
codes the period of the cycle which is also given by the framed numbers for trajectories 
falling in rhombuses. The circled numbers stand for the periods of the cycles falling 
on dashed lines. 
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Figure 8. (Color online) Bifurcation diagram (bottom panel) and Lyapunov exponent 
(upper panel) obtained using b = 2.01 and r = 4. This plot presents a zoom of Figure 
[71 The inset in the upper panel allows to appreciate the plateau in A. 
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Figure 9. (Color online) Bifurcation diagram (bottom panel) and Lyapunov exponent 
(upper panel) obtained using b = 2.01 and r = 4. This plot presents a zoom of Figure 
[8] in the region located by the arrow. It exhibits interior crises. 




Figure 10. (Color online) Bifurcation diagram (bottom panel) and Lyapunov 
exponent (upper panel) obtained using b = 2.01 and r = 4. This plot presents a 
zoom of Figure m The inset in the upper panel allows to appreciate the plateau in A. 
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Figure 11. (Color online) Detail of the first 350 Harter curves obtained using b = 2.01 
and r = 4. Dashed (red) lines between crossings of harter curves correspond to unstable 
fixed points. 
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Figure 12. (Color online) Left panel: Lyapunov exponent as a function of c for three 
values of r. Right panel: Collapse of the same three curves as a function of ^ 
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Figure 13. (Color online) Bottom panel: Bifurcation diagram as a function of c, with 
b — 2.1 and r — 4. The arrow locates the narrow window zoomed in the upper panel 
where A is also shown. 




Figure 14. (Color online) Bifurcation diagram and Lyapunov exponent as a function 
of c. This particular plot was obtained using b = 2.2 and r — 2. The system is not 
close to the critical point 6 = 2. 
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Figure 15. (Color online) Bifurcation diagram and Lyapunov exponent as a function 
of c, with b = 3 and r = 2. 




Figure 16. (Color online) Crisis-induced intermittency, with r = 4 and b = 2.2. Two 
trajectories from both sides of the critical value c* = 0.88411175... are compared. 
Solid dots (red) correspond to a trajectory with c — 0.8845, where the attractor splits 
in narrow bands. Empty squares (black) correspond to c = 0.8840, namely in the wider 
region of the attractor. The inset is a zoom of the points close to x = 2.5 and allows 
to illustrate the non-regular character of these trajectories. 
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Figure 17. (Color online) Characteristic length r versus |c — c*| in log- log scales. 
Each dot was obtained from 100 initial conditions and series of 10^ iterates, with 
c* = 0.88411175, b = 2.2 and r = 4. The slope of the linear fit is 7 = -0.89 ± 0.02. 



